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Abstract. 


A digital computer program CUBIC has been developed 
for on-line differentiation of analog voltage signals. The 
program accepts voltage records of a time function and yields 
its derivative after one program cycle time of 4.2 msecs. The 
velocity is obtained by employing a least mean squares cubic 
fit technique. 

The routine is intended for experimental work either 
as a data reduction tool or as a control signal for a closed 
loop experiment. The program can be implemented on a PDP-8 
digital computer with one digital to analog converter channel 
and one analog to digital converter channel. 








i. 


The derivative (hereafter called the velocity) of 
an analog signal or a continuous time function (hereafter called 
the position) is of interest in many classes of research. 

This information can be obtained by a variety of techniques. 

The approach used in this program is to fit a cubic curve 
through the last 33 samples of data using a least mean squares 
fit technique and taking the slope of this curve at the latest 
data point. The program developed produces the least noise of 
all differentiation techniques considered (a six point difference 
method, pseudo differentiation with a transfer function as/(s+a), 
and a parabolic least mean squares technique) and yet is applicable 
over a wide range of frequencies: d.c.'to 8 Hz. The cubic fit 
algorithm, together with examples of its operation, is described in 
reference 1. 


Operation of t he program . 

The program is designed for operation in a hybrid 
installation. The hybrid equipment used consisted of a Digital 
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1.5 usee, cycle time with Extended Arithmetic Element Type 182) , 
connected to a GPS-290T analog computer with ±10 volts reference 
voltage. An analog clock generated the sampling signal for the 
A/D and D/A conversions at a rate of 160 per sec (or every 


r 


6.25 msecs .) 


Storage. 

The digital progreim is compiled as an eight page sub¬ 
routine and is listed in Appendix A. The program occupies locations 
(20 ,34), (47,111), 126, 127, (5233,5334 ) and (5600,7567). The auto¬ 
index 17 is used and the 33 data points are stored in locations 
(6G00 f 6040). Double and single precision multiplication routines 
for signed numbers are required by the routine (ref.2). 
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Precalculated constants . 

The digital program has been compiled assuming 33 

data points. Assume N is the number of data points. The 

3 

constants N, (N-l)*(N-2), and (N-l) are placed in locations 

2 

6377, 6376 and 6452 respectively. The constant (N-l) is formed 
in instructions 6311 and 6312. Choices other than 33 data 
points are permissable. 2 r +l, where r is an integer,is an 
advantageous choice given the binary basis of digital computers. 

The values of recursion formulas used in the program, 
for a unit step in position of duration equivalent to the N data 
points, are stored in locations (7534, 7540). These values are; 
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N-l 


N-l 
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i=0 

N-l 

E 

i=0 

(~i) 

respectively 


. v 2 


N-l 


. x 2 


The upper part of E (-i) is formed at location 7456, and 

i=0 

N-l N-l 

the upper part of £ —(i—1) at location 7476. E (2i-l) 

i=0 i=0 


is computed by instructions (7520, 7524). 

The velocity estimation equation prior to digital 
filtering may be written as 

aj = XVEL = [ lOct-jA 11 + 10a 2 B n /At + 10a 3 C n /At 2 + 10a 4 E n /At 3 ] * 


[ l/10At ] (1). 

4 

At is the sampling interval. The constants 10a^, 10a ? , 10a 3 (2 ), 
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and 10a 4 (2 ) are precalculated in double precision arithmetic and 
located at (6750, 6757) . The constants a ± (i=l,4) and the recursion 
formulas A , B , etc. are defined in reference 1. Both parts of 
these double precision fractions must be expressed as a signed 
fraction. For example, the fraction 
0.12515777 g 
is expressed as 

0.1252g and -0.00002001g 

and stored as 

0.1252 and -0.2001 . 

The data points may be stored in any N consecutive 
core locations. Locations 7565 and 7566 must contain the last 
and first data storage locations respectively. The order of 
storage is inverted with respect to sampled time. 

Calling sequence . 

The program is called as a subroutine. A typical 
calling sequence including the A/D and D/A conversions is listed 
in Appendix A. X2 is the filtered estimate of the velocity. 

Flow of operations . 

Figure 1 is a flow diagram describing the sequence 
of operations. 

the AC must be clear at entry, the link status is 
irrelevant. 

— to restart estimation, START (location 52) must 
contain a negative value. 

the jump is made to the subroutine with the latest 
data value at MTHETA (location 55). 

— the recursion formulas A n , B n ,.etc. are updated. 
a l^t is calculated in triple precision arithmetic, 
the unfiltered velocity estimate, a^ , is calculated 
by shift operations, simultaneously a scaling by 2 
can occur. 

the jump is made to the subroutine LIMIT. 
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— the jump is made twice to the subroutine FILTER 

— the jump is made to the subroutine LLAGS 

— the filtered velocity estimate is stored at X2 
(location 126) 

— the jump is made back to the main program. 

— the AC is clear at exit, the link status is 
undetermined. 


The calculation of a^ from a^At . 


In the flow diagram (figure 1) B n /At is written as 
B(n) , similarly for C(n) , etc. At location 6732 the first 
factor within square parantheses in equation 1 above is con¬ 
tained in four locations designated FVH, FVM, FVL, and FVLL. 

The binary point is between bit 1 and bit 2 of the register FVM. 
(The correspondence between the analog voltage reference levels 
and the digital signed fraction system is 0.3777 o = 10 volts on 

O 

the above cited hybrid installation.) 

lOAta^ is independent of At, thus the program operator 
may choose any scaling factor or sampling rate to obtain the 
velocity estimate in operations subsequent to that of location 
7004. Multiplication by l/10At is affected by choice of the 
variable SCALE (location 7045) . The subroutine SHIFT shifts the 
contents of FVH, FVM, FVL, and FVLL left or right to obtain the 
required multiplication. In the present application one shift 
left is required to place the decimal point between bit 0 and 1 
of register FVM and as At equals 6.25 msecs. 4 further shifts 
are required for multiplication by l/10At. Thus the total 
content of SCALE is +5. It is recommended that sampling rates 
which are a binary multiple or divisor of 160 per sec. be used. 
Storage of the velocity estimate on return from the subroutine 
SHIFT is at location 7324 (ULOAD). If the estimate exceeds 
±10 volts the program halts at location 7016 or 7020. 


The subroutine FILTER . 

Differentiation routines are very sensitive to high 
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frequency noise. To reduce this noise two equal cascaded simple lag 
filters are used. The input to the filter cascade is ULOAD and the 
output is XVEL (location 56). FILTER is one of the lag filters, 
(a/(s+a)), with a cut-off frequency determined by the subroutine 
LIMIT. The digital filter is based on the state space representation 
of a simple lag: 


uft^ - $ ( t i/ t -j L _i) u (t i _,) + a f 1 $ (t • ,x) c( x) dx, 

i-l 1 

where u(t ± ) is the output of the filter and c(x) is the input. 
As a difference equation 


t. 

u(i) = $ (At) u (i-l) + ac (i-l) / 1 $(t.,x)dx 

t. i-l 1 

The terns $(At) and a / 1 $(t. ,x)dx are stored at PHIX and 

i-l 1 

PHIV respectively (73 76 and 7375) . The program is compiled with 
a preset at 20.4 Hz. 


The subroutine LIMIT 

A flow chart for the subroutine LIMIT is shown in 
figure 2. The function of the subroutine is to determine an 
appropriate lag constant a for the subroutine FILTER. This 
is accomplished by determining the maximum velocity (i.e. of 
XVEL) and maximum input position (i.e. of MTHETA) to obtain 
an estimate of the fundamental input frequency. Based on this 
estimate the lag constant is set as follows: 

8 when oxl, 

8ti) when 64Xo>l, 

512 when go>64 

where u) is the frequency in radians per second. 

If the velocity estimate was a single sinusoidal time 
function with no attenuation with respect to the true velocity 
of the input position sinusoidal time function, the a chosen by 
the subroutine LIMIT would introduce a phase lag of 11.4 degrees 
into the velocity estimate. In practice the velocity estimate is 
noisy; the phase lag introduced is about 10 degrees in the range 
u<64. A fixed lag constant may be used by setting the variable 

HOLD (location 6161) to a/8 and noting that the binary point lies 
after bit 11. 
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Bode plots . 

Bode magnitude and phase plots for CUBIC are shown in 
figures 3 and 4. The -3 db. point, with respect to a pure 
differentiator, is at 9 Hz., whereas the 45 degrees phase point is 
at 4.5 Hz. X2 is the filtered estimate of the velocity, being 
the output of the subroutine LLAGS. LLAGS is a lag-lead digital 
filter which corrects for high frequency amplification in the 
range 2.5 to 9 Hz. This amplification is a direct result of the 
least mean squares cubic fit procedure. 
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MTHETA is a new data point. 


6200 C STARTUP ) 


6210 

B (n) B (n+1) 

6235 

D(n-l) «- D(n) 

6267 

C (n) «- C (n+1) 

6317 

-G(n-l) 4- -G(n) 

6366 

E(n) + E(n+1) 

6603 

A(n) A(n+1) 


7004 


6453 


a 1 (n) = 7324 



7421 


All data buffer «- MTHETA 
If MTHETA = V 
A(n) 41V 
B (n) «- -1020V 
C(n) •+■ 26260V 
D(n-l) 2000V 
-G(n-l) «- -25240V 
E(n) «- -1040400V 
«- 0 



Figure 1. Flow Chart for CUBIC. 
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Flow Chart for Subroutine LIMIT 

























































Fig. 2 . Bode magnitude plot for CUBIC . 
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